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A simple new algorithm for the calculation of two-particle 
Bose-Einstein correlations from classical event generators is 
derived and discussed. 

PACS numbers: 25.75.+r, 07.60. ly, 52.60.+h 

: 

0*\ For ultrarelativistic heavy- ion collisions, two-particle 
T3ose-Einstein correlations between identical pions or 
ions provide a unique possibility to reconstruct the ge- 
ntry (size, temporal extension) and dynamics (col- 
sctive expansion flow) of the source at the point of 
r4*adron freeze-out. (For a recent theoretical review see 
) This reconstruction is, however, not completely 
model independent; it requires the use of "reasonable" 
s§purce parametrizations (e.g. [0-^|) whose parameters are 
then fixed by a simultaneous analysis of single-particle 
T^nomentum spectra and two-particle momentum correla- 
Oions g|. 

^J" Invaluable help for the selection of "reasonable" source 
i^parametrization comes from microscopic event generators 
^e.g. VENUS ||, RQMD § or ARC §) which gener- 
~*2i£e the phase-space distribution of hadrons at freeze-out 
+2!om a dynamical Monte Carlo simulation of the (clas- 
^cal) kinetic phase-space evolution of the collision zone. 
^Unfortunately, it was recently pointed out |||| that the 
mrect computation of two-particle correlation functions 
n3om classical kinetic codes ^j-Q is fraught with severe 
m Conceptual and practical problems. These can be sim- 
ply understood by starting from the general relation ]l0| 
thjetween the 2-particle correlation function C(q, K) and 
5$ie (real) "emission function" (one-particle Wigner den- 
sity at freeze-out) S(x,K), 
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(p a ,E a ) and (pb,Eh) are the 4-momenta of the two ob- 
served particles. Eq. ([la]) neglects final state interactions 
which we will leave out in this Letter in order to concen- 
trate on the principal issues. 

The numerator in the second term of Eq. ([la|) can be 
rewritten as 



d 4 xd 4 yS(x,K)S{y,K) 
x cos(q-(x - y)) . (2) 



The problem is the construction of the Wigner density 
S(x, K) from the output of the event generators. The lat- 
ter consists of a set of phase-space points (xi,pi) denoting 
the (on-shell) momenta pi and points of last interaction 
Xi of the produced particles. According to ( |lb| ) K is the 
average of two on-shell momenta, but not itself on-shell, 
K a \lm 2 + K 2 . Therefore S(x, K) cannot be directly 
related to the phase-space density generated by the dis- 
tribution of points (xj,pj). To overcome this difficulty 
one ususally imposes j[l],[l2) the "smoothness assump- 
tion" S(x, Pa + Pb ) w S(x,p a ) w S(x,p b ) and rewrites the 
expression (|2|) as 

Num(g,if) = J d 4 xd 4 y S(x,p a ) S{y,p b ) 



x cos((p a -p b )-(x - y)) 



(3) 



One now identifies S(x,p) with the classical output dis- 
tribution from the event generator, 



N 



Sclass (X,p) = y^^ (4) (x - Xj) 8^>{p- pi) . 



(4) 



where N is the total number of pions of a given charge 
in the event. In the widely used Pratt algorithm fll2|] 
the expression which results after inserting (|j) into (|3|) is 
simulated by the ad hoc prescription 

Num(q,K) i — ► ^ C0S ((Pi ~ PjY( x i ~ x j)) ■ ( 5 ) 

Here "bin" denotes a small bin in (q, K) withp^ — pj w q 
and (pi +pj)/2 w K. (In practice the bin size depends 
on event statistics.) 

The prescription (Q) has two severe problems: first, the 
positivity of Num(g, K) got lost between Eqs. (Q) and 
([}]) when making the "smoothness assumption" . It was 
pointed out in Refs. |l3|,[l4| and practically demonstrated 
in Ref. |)j that for sources with strong x-p-correlations 
(e.g. rapidly expanding sources) this can lead to un- 
physical oscillations of the simulated correlation function 
around unity. Second, the intuitive substitution law (||) 
is formally incorrect and results in a wrong selection of 
contributing pairs as well as an incorrect weighting factor 
for each pair. 

To prove this last point let us write in (^|) 



S(x,K) S(y,K) = J d 4 P 1 d 4 P 2 S(x,P 1 )S(y,P 2 ) 
x S {4 \P 1 -P 2 )S^ (l 



P1 + P2 



(6) 
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Inserting the classical expression (J4j) one finds instead of 
® 
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Here "bin (if, e)" denotes a small bin around if with 
width e in each of the four directions. The prescription 
( fflb| ) can be derived rigorously and directly by first gen- 
crating from Eq. (^) a piecewise constant function ( "his- 
togram" ) through "binning" , 

S(X,K)= / d^pSdaesiXtP) 
Jbin(K,e) 
rK+t/2 

d 4 pS c i ass (x,p) , (8) 

'K-e/2 

and then inserting S(x, if) into Eq. (|J). Eq. @ is a tech- 
nical step required by finite event statistics; in practice, 
e should be chosen as small as technically possible. Note 
that the selection of pairs in J7b| ) differs from the one 
in (||): for given if the algorithm ([?]) selects pairs with 
pi w w if, independent of the value q at which the 
correlation is to be evaluated. For different values of q at 
fixed if, the correlator is obtained by weighting the same 
set of pairs with different weight factors cos( q-(xi — Xj)) 
which depend only on the spatial coordinates, but not 
on the momenta of the particles in the pair. This is con- 
sistent with the expectation from Eqs. (|la|,||) that the 
measured ^-dependence of the correlator gives access to 



the distribution of relative distances X; 



in the source 



(at fixed if). Since the steps from Eq. (|J) to Eq. (0) 
involve only identical transformations (the difference be- 
tween ( |7a| ) and (7b) arising only from a different choice 
of emission functions (^) resp. ([§)), they preserve the 
positivity of the second term in (|la). Please note that in 
contrast to (||), (0) is a continuous function of q, i.e. no 
binning in q is required. 

We have checked for the simple analytically solvable 
model presented in || that the algorithm (]?]) indeed al- 
lows to reconstruct from a classical Monte Carlo simu- 
lation the correct analytical expression for C(q,K); in 
particular it removes the unphysical oscillations found in 
||. In Fig. 1 we show various approximations for the cor- 
relation function for a 1-dimensional source in z-direction 
with emission function 



S(z, t; if) = e~ z2 / R2 S(K - az) S(t) , 



(9) 



with i?=10 fm and a— 0.02 GeV/fm. This classical source 
features perfect z-if-correlations, K=az, which are, of 
course, quantum mechanically forbidden (see below). For 



the source (^|), the exact correlator (1) can be calculated 
analytically, yielding C(q) = 1 + exp[g 2 /(2a 2 i? 2 )] (solid 
line in Fig. 1). The pathological rise of C(q) above the 
value 2 at q=0 is due to the violation of the uncertainty 
relation between z and p z by the model @; we selected 
this model because we believed that such a feature may 
be particularly difficult to reproduce in an event genera- 
tor. Indeed, reconstructing the correlator from a Monte 
Carlo simulation of (|9|) via the Pratt prescription (||) 
yields the long-dashed line in Fig. 1 ||; it can be calcu- 
lated analytically from Eq. (||) as C(q) = 1 + cos(<7 2 /a). 
This is always less than 2, but oscillates wildly around 1, 
becoming even at regular g-intervals. This contradicts 
the positivity of the second term in (1) and is due to 
the smoothness approximation (||). - The two remain- 
ing lines in Fig. 1 show results from the same Monte 
Carlo simulation of (0) but reconstructing the correlator 
through the new algorithm (7b) instead of (0). For the 
dot-dashed line the bin width e was chosen as e = 10 
MeV, for the dotted line as e = 5 MeV. In both cases 
the simulated result deviates from the exact one (solid 
line); this is not a failure of the algorithm, but a result 
of the binning procedure (||) applied to the source (||) - 
a purely technical step required by finite event statistics. 
As seen, the discrepancy decreases with decreasing bin 
width e. 
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FIG. 1. Two-particle correlation function for the model 
source (^). Different curves show different approximations as 
described in the text. The statistical errors of the simulation 
are below the line widths. 

While Eq. (^) thus solves the technical problems of 
the prescription (||), it does not address the principal 
physical problem that the classical distribution (^) is 
not a valid Wigner density since it violates the uncer- 
tainty relation by simultaneously fixing the coordinates 
Xi and momenta pi of the emitted particles. It has been 
repeatedly suggested [|l5] that this can be remedied 
by replacing the sharply localized (^-functions in (^) by 
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minimum-uncertainty wave packets. In the remainder 
of this Letter we will discuss how such a procedure will 
modify the algorithm for calculating single-particle spec- 
tra and two-particle correlations from event generators, 
thereby rendering it quantum mechanically consistent. 

Let us start from the folding relation for the emission 
function, derived in Ref. |l3j within the covariant current 
formalism p8| : 

S(x,K) = J d A zd i Qp{x- z,Q)S a {z,K - Q) . (10) 

Here 

S (x,p) = J d^ve-^j* (x+ f) jo (x- |) (11) 

is the Wigner density associated with an elementary 
source current amplitude jo(x), taken below as a Gaus- 
sian wavepacket, and p(x,p) is a classical phase-space 
distribution for the centers Xi and (average) momenta pi 
of these wave packets, here taken as 



N 



p(x,p) = J2t i4 Hx~x l )6^(p-p l ). 



(12) 



i=i 



For the elementary source amplitude we make the ansatz 



(13) 



this source emits a Gaussian wave-packet with width pa- 
rameter a at "freeze-out time" a; . With this ansatz the 
elementary Wigner density Sq becomes 



x 



So(x,p)=8(7rcT 2 )t\M\ 2 5(x {, )exvl [ --2-a'p 2 j . (14) 
Inserting this and ( [L2] ) into (10) one finds 

S(x,K)=Af>jr S(x°-x°) exp(- (3; 'f )2 ) 
i=l \ a / 

x e X p(-a 2 {K- Pi ) 2 ) , (15) 

with Af' = 8(na 2 )i\J\f\ 2 . This generalizes the classi- 
cal ansatz (^) into a quantum mechanically consistent 
source Wigner density; with the free parameter a one 
can choose the relative degree of localization in coordi- 
nate space (cr — > 0) or momentum space (a — ► oo), always 
preserving Ax • Ap = h. 

The single particle spectrum, which occurs in the de- 
nominator of the second term of (P~g), is now given by 



E 



dN 



d 3 Pa 

Af" = (2tr7 2 ) 2 |/v1 2 , 

Vi(p a ) = exp(-<T 2 (p a - Pi) 2 



N 

d i xS(x 1 p a )=AT"Y,v l (Pa), (16a) 

i=l 



(16b) 
(16c) 



and similarly for pb. It is normalized to the total num- 
ber N of pions in the event; this fixes the normalization 
constant Af above. The exchange term (0) in the two- 



particle spectrum (with q and K defined in Eq. (lb)) is 
similarly derived as 



Num(q, K) = {Af") 2 cxp(- \a 2 q 2 



N 

(17a) 



j(q,K) = Vi(K)vj(K) cos(q-(xi - Xj)). 



(17b) 



The normalization Af" drops out in the correlator (fj"a|). 
Eq. dl7) should be compared with the classical expres- 
sions (|) and (0). Like (Q) (but contrary to (||)) it is posi- 
tive definite and thus free of spurious oscillations around 

0. The sum in ( |l7a|) is now over all pairs the sharp 
restriction to the bin "bin(i^, e)" in (0) is replaced by the 
Gaussian weight factors Vi{K) and exp(— \a 2 q 2 ). Please 
note, ho weve r, that the former also occur in the new def- 
inition ( |l6a|) for the single particle spectrum, and must 
be kept in both places for consistency. Eq. ( 17b ) shares 
with Eq. (0) the property (which was already discussed) 
that the cosine weight factor for each pair (i,j) depends 
only on xi — Xj, but not on pi and pj. By combining 
Eqs. ( filf ) and ([17]) it can be shown that now the corre- 
lator C(q, K) is always between 1 and 2, i.e. that the 
pathological rise above 2 shown in Fig. 1 cannot happen 
for an emission function which respects the uncertainty 
relation. In Ref. |L9| the results @ and ([j]]) (but not 
(^|)) were derived with different methods, including finite 
multiplicity corrections. 

Expressions (|l^) and (jlj]) depend on one free param- 
eter, the Gaussian width a of the wave-packets. It is 
instructive to discuss the two obvious limits, a —>■ 
and a — > oo. For a = the elementary sources are 
sharply localized in space (cf. Eq. (|l^)); as a result, the 
single-particle momentum spectrum (|l6|) is completely 
flat. Furthermore, one sees from Eqs. ( |l7|) that the corre- 
lation function becomes if-independent, even for expand- 
ing sources of type (|l2|) where the Xi and pi are strongly 
correlated. Both features are clearly unrealistic. In the 
opposite limit, a — > oo, the Gaussian smearing facto rs in 
the single particle spectrum ( [l6| ) disappear, and Eq. (16a) 
degenerates to a sum over ^-functions; this is the usually 
employed algorithm for computing single-particle spec- 
tra from event generators with classical particle propa- 
gation The two-particle correlation term (|l7|), on 
the other hand, is then sharply concentrated at q = 0, 

1. e. the correlator C(q, K) drops from 2 to 1 over a q- 
range of order 1/er. This translates into a source radius 
~ a and reflects the diverging spatial extension of the 
elementary wavepackets in this limit, irrespective of the 
(localized) spatial distribution p{x,p) of their centers. 

It is thus clear that in practice tr must be kept fi- 
nite, but non-zero. As pointed out in Jl5|-JT7t this im- 
plies a broadening of the single-particle momentum dis- 
tributions relative to the one derived from the classical 
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phase-space distribution p(x,p) of freeze-out points. Us- 
ing the algorithms (16 1?]) for a quantum mechanically 
consistent computation of spectra and Bose-Einstein cor- 
relations from classical event generators thus requires a 
retuning of the codes to elementary e + e~ and pp colli- 
sions, using the same algorithms there. 

This last step imposes rather restrictive limits for the 
value of a Jl5[ . Since the pions from high energy pp 
and e + e~ collisions have an average transverse momen- 
tum (p±_) ~ 0.35 GeV, the Gaussian width in momen- 
tum space of the elementary wave packets must be be- 
low this value. This implies a > 0.5 fm. On the other 
hand, the effective source radii for such collisions ex- 
tracted from 2-particle correlations are of the order of 
only 0.8 — 2 fm |pofl . This implies er < 1 fm. It was 
already pointed out in the pioneering GGLP paper |22j 
that this value roughly agrees with the pion's Compton 
wavelength. This suggests the following interpretation 
of the measured pion spectra and Bose-Einstein corre- 
lations from elementary hadron-hadron and e + e _ colli- 
sions: each elementary collision produces a small number 
of elementary Gaussian wave packets with width a, 0.5 
fm < fj < 1 fm. For a ~ 0.5 fm the width of the Gaus- 
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sian momentum distribution in (16c) nearly exhausts the 
measured the measured single-pion spectra thus 

reflect mostly the intrinsic momentum distribution of 
the elementary pion wave packet. On the other hand, 
the HBT radius i?™ T = a /V2 corresponding to the q- 



Gaussian in (17a) nearly exhausts the values extracted 
from two-particle correlation measurements; the latter 
thus mirror the intrinsic spatial width of the wave pack- 
ets. There is very little additional room in the data for 
random (thermal) motion of the wave packets relative to 
each other, nor for their spatial distribution over a larger 
volume. If the total source is much bigger than 1 fm, 
it must expand very rapidly with homogeneity regions 
which are not much larger than the size of an elementary 
wave packet. In this sense, pion spectra from e + e" and 
pp collision measure the smallest sources compatible with 
the uncertainty relation. 

In summary, we have derived a new algorithm for the 
computation of single-particle spectra and two-particle 
correlations from classical event generators. The classi- 
cal algorithm (^) removes the recently discovered defi- 
ciencies of the presently employed Pratt algorithm (|B|). 
The quantum mechanical algorithm ( |u^ , |l7| ) additionally 
ensures that the uncertainty relation is not violated. We 
also showed how the free parameter a in the latter al- 
gorithm can be fixed from elementary e + e~ and hadron- 
hadron collisions. 
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